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Abstract. 

We present a method to numerically compute accurate tunnelling rates for a Bose- 
Einstein condensate which is described by the nonlinear Gross-Pitaevskii equation. 
Our method is based on a sophisticated real-time integration of the complex-scaled 
Gross-Pitaevskii equation, and it is capable of finding the stationary eigenvalues for 
the Wannier-Stark problem. We show that even weak nonlinearities have significant 
effects in the vicinity of very sensitive resonant tunnelling peaks, which occur in the 
rates as a function of the Stark field amplitude. The mean-field interaction induces a 
broadening and a shift of the peaks, and the latter is explained by analytic perturbation 
theory. 
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1. Introduction 

Quantum dynamics often is intriguing and counter-intuitive. A prominent example 
thereof is the locahsation of a wave packet in a spatially periodic lattice induced by an 
additional static force: the force can turn an extended Bloch wave (which is a solution of 
the Schrodinger equation with a periodic potential fX]) to a wave packet which oscillates 
periodically in (momentum) space p. While conceptionally simple, this well- know 
Wannier-Stark problem is complicated from the mathematical point of view because 
the system is open, i.e., unbounded, and any initially prepared state will, in the course 
of time evolution, decay via tunnelling out of the periodic potential wells |21 E] • 

Starting from the Bloch bands of the unperturbed problem (i.e., without the static 
field F = 0), the decay can be attributed to tunnelling from the ground state band to the 
first excited energy band. The celebrated Landau-Zener theory predicts an exponential 
decay rate (see, for instance, [HE] for introductory reviews): 

r(F) oc Fe-^ , (1) 
§ Corresponding author's e-mail: saw@df.unipi.it 
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where h is proportional to the square of the energy gap between the two lowest energy 
bands. For experiments with cold atoms, i.e., the scenario on which we focus in this 
paper, the wave packet decays very quickly by successive tunnelling events, once it 
has tunnelled across the first band gap. This is due to the much smaller gaps of the 
higher energy bands in a sinusoidal potential [HI Ej. The Landau-Zener formula (|T|) 
cannot account for the interaction of the Wannier-Stark levels at adjacent potential 
wells. Between such adjacent lattice sites nearly degenerate Wannier-Stark levels repel 
each other, which leads to a strong enhancement of the tunnelling decay. These resonant 
tunnelling events result in pronounced peaks in the rates as a function of the inverse 
field amplitude on top of the of the global exponential decay described by ((H) [51 IH]. 

Figure ^a) shows two Wannier-Stark levels on each lattice site. The levels within 
either of the two ladders are separated by mFdi in energy, where di = 27r denotes the 
lattice period and the integer m counts the number of sites in-between two energy levels 
of the same ladder jH |2] . The decay rates for non- interacting particles in the periodic 
potential V{x) = Vq sin^(a;/2) +Fx can be computed from the Wannier-Stark spectrum, 
e.g., by using the numerical method described in [S]. Figure [Tfb) presents the rate F as 
a function of 1/F. The maxima occur when mFdi is close to the difference in energy 
(El — Eq) between the first two energy bands (averaged over the fundamental Brillouin 
zone in momentum space) of the unperturbed {F = 0) problem [21 E]- The actual peak 
positions are slightly shifted with respect to this simplified estimate [marked by arrows 
in Fig. ^b)], owing to a field-induced level shift close to the avoided crossings of the 
levels P^. 

Exceptional experimental control is possible nowadays with Bose-Einstein 
condensates (BEC) whose initial conditions in coordinate and momentum space can 
be adjusted with unprecedented precision. With the help of a BEC, sensitive tunnelling 
phenomena were studied in time-dependent systems ^Uj, as well as in static potentials 
[TT] . Here we are interested in tunnelling in the Wannier-Stark problem where 
the impact of the intrinsic atom-atom interactions in the BEC has been studied in 
several recent experiments [3 El El E]- In those experiments, the difficulty in 
understanding quantum transport processes, such as coherent tunnelling, originates 
from the complex interplay between classical transport in the underlying phase space, 
quantum interference effects, and the many-particle interactions. 

In a typical experiment with a BEC, where the number of atoms in the condensate 
is large and the atom-atom interaction is rather small, the Gross-Pitaevskii equation 
(GPE) describes the condensate in very good approximation ^21- Recently, some of 
us proposed a concrete experimental scenario to measure the impact of a mean-field 
interaction potential (which in the GPE takes into account of the atomic collisions) on 
the tunnelling in the Wannier-Stark problem jTHj. More specifically, the interaction- 
induced modification of the resonant tunnelling peaks was studied, and it was found 
that the peaks [such as the ones in Fig. Wljo)] are washed out for a large enough - 
but still experimentally feasible - interaction strength. As discussed in jTH], for a finite 
mean-field nonlinearity, the concept of decay rates is not as well defined as in the case 
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Figure 1. (a) Schematic skctcli of nearly degenerate Wannier-Stark levels (thin line: 
ground state levels; thick line: first excited levels in each well) in a potential of the 
form V{x) — Vq sin^(x/2) + Fx. (b) tunnelling rates F for Vo = 2 as a function of the 
inverse Stark field amplitude 1/F. The resonant tunnelling leads to the pronounced 
peaks, which lie approximately at F « {Ei — i?o)/(27rm) (with integer m). These 
estimates (marked by the arrows) are slightly modified by field-induced level shifts. 

of non-interacting particles. The reason is that the weight of the nonlinear term, which 
is proportional to the condensate density, varies in time. Hence, the probabihty that 
an initially prepared state will stay in the preparation region does not follow a simple 
exponential law. In other words, the nonlinear interaction decreases as the condensate 
escapes via tunnelling and, as a consequence, the decay rates can be defined only locally 
in time. 

In this paper we want to discuss how the problem of defining proper decay rates can 
be solved. One way to define a time- independent, global tunnelling rate is to renormalise 
the density of the condensate in the preparation region (e.g., in the central potential well 
of the periodic lattice) continuously, such that the average density remains constant in 
time. For experimentally realisable nonlinearities El El E| , this approach results in 
a mono-exponential decay of the survival probability, and in consequence in a reasonable 
definition of the decay rate. The corresponding resonance states are characterised by 
their stationary asymptotics, much in the same way as the stationary solutions of a linear 
scattering problem (i.e, described by the linear Schrodinger equation) 0. We solve 
numerically the well-posed problem of finding the resonance states, using the method 
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of complex scaling. The theoretical background to treat the nonlinear interaction in 
the GPE in a consistent way was recently laid in ^H]- We use a modified version of 
this method, with crucial extensions on the algorithmic side, which proved necessary to 
stabilise the computations for more complicated potentials than the single- well potential 
exemplarily treated in |18j . 

We review the defining equations of the nonlinear Wannier-Stark problem and the 
complex scaling technique of ^H] in the following section |21 where we also describe our 
numerical algorithm in detail, c.f., subsection 12 .21 Section IHl presents our central results 
on the decay rates in the vicinity of the resonant tunnelling peaks for experimentally 
relevant nonlinearities. Section |3] finally concludes the paper. 



2. The nonlinear Wannier-Stark problem 

We use the one- dimensional GPE to model the temporal evolution of a BEG loaded into 
a spatially periodic optical lattice potential and subjected to an additional static force 
F: 

1 (9^ _ . o /x^ 



+ K)Sin2 +i7^a; + ^|^(^^t)|2 



^{x,t). (2) 



2(9x2 

ilj{x,t) represents the condensate wave function, and we used the dimensionless 
quantities Vq = Vsi/Eb,F = Fa^idL/ {2tiE^), g = gsid^N/ {27!- E^). The characteristic 
length scale is the lattice period rf^, i.e., x = xsi^tt/cIl, the Bloch energy is Eb = 
(irh/dL)'^ /M, with atomic mass M, the number of atoms N, and the (from three to one 
spatial dimensions) rescaled nonlinearity parameter ^fgi (see [12] for a definition of gsi, 
where also the regime of validity of the one-dimensional approximation is discussed in 
detail) . 

Since V{x) = Vosin^(x/2) + Fx — oo for x — oo, any state initially prepared 
in the optical lattice will escape via tunnelling. We search for the resonance state ipg 
which solves the stationary version of Eq. (j2|) 

H[i^9]i^g = Egij, , (3) 
for the eigenvalue Eg = fig — iTg/2, and the Hamiltonian 

m = -l^ + V{x)+gmx)\' (4) 

To render the problem posed by Eq. (jH)) meaningful we demand that the condensate 
wave function remains normalised around the initially prepared state, i.e., around x ~ 0: 

dx \iljg{x)f = 1 . (5) 

The boundaries x„ must be chosen in a reasonable way, and we chose x„ = tt (so the 
probability to stay in the central well around x = remains one ^^). We verified that 
slightly different choices of the boundary vr < x„ < 3tt/2 led to eigenvalues which did 
not change on the significant digits given below in section |21 
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We discuss now the renormalisation condition (0) and its consequences. In practice 
such a condition may be reahsed by the presence of a source term which constantly 
supphes a condensate flow j25- Experimentally such a scenario could be achieved by 
constantly reloading the central well with coherent BEC matter. Transport experiments 
of such kind could be realised with the help of optical tweezers [22], atomic conveyer 
belts [221 5 or microscopic guides for ultracold atoms [211. realisation may introduce 
additional modifications in the temporal evolution of the decaying system, which go 
beyond our simplified assumption of renormalisation. Such modifications, e.g., the 
relaxation of added particles in the periodic lattice potential, depend on the specific 
realisation. We expect, however, that the asymptotic decay will be hardly affected by 
such processes, the time scales of which should be relatively short and of the order of 
the period of oscillations in the potential wells. 

On the other hand, if the condensate wave function tunnels out of the central well 
without sudden changes of its shape, the time-dependent atomic population N{t) inside 
the well decays according to the relation [THl I2S] 

'-^ - -T,,m) . (6) 

Assuming that the decay rate adiabatically adjusts itself to the time- dependent value 
Tg{t), with g{t) oc N{t), Eq. © can be solved for a given initial number of atoms A^(0) 
in the condensate. Knowing the "local" rates Tg(^t) for Q <\g\< 1(7(0)1 allows us then to 
compute the actual survival probability in the central well, which in [T^ was obtained 
differently by a brute force integration of the time-dependent GPE ([21). 

We emphasise that the setup studied in Ref. [12] bears some crucial differences 
to the problem posed here, which is based on condition In ^B] the short time 
behaviour of the relaxed ground state (for F = in the periodic potential and in the 
presence of additional harmonic confinements) was predicted for the three dimensional 
Wannier-Stark problem. The approach presented here is capable to determine, via 
Eq. ([nj, the decay only for single resonance states according to the above arguments. 
Although such resonance states are typically distributed over many lattice sites, they do 
not provide a prediction for the decay of a general initial state (which could be composed 
of contributions from many adjacent wells), simply because the superposition principle 
does not apply for the nonlinear GPE ([2)). 

In this paper we want to compute directly the precise decay rates Vg of a single 
resonance state using the complex scaling method, which is described in the following 
subsection. 

2.1. Complex scaling 

For the linear problem with (7 = 0, one of the standard techniques to compute resonance 
states numerically is the complex scaling method (which goes back to and is 

reviewed, for instance, in [2Z])- Applying the renormalisation condition ^ allows us to 
use this method to find the stationary eigenstates and the corresponding eigenvalues. 
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see Eq. 0. Without this condition, the nonhnear interaction term would vary in time, 
and a stationary state would not exist because of the tunnelling decay. 

An additional problem when dealing with the nonlinear term in the GPE arises from 
the method of complex scaling itself. The problem of defining the complex conjugate 
of the wave function ifji^x) is described in [IHIESI, and was solved in [THj. Usually, the 
scaling transformation is defined as follows 

^(x) ^ /(x) = R{e)^{x) = e'^l^^[xe'^) , (7) 

where the pre-factor is just a phase depending on the dimensionality of the problem (here 
we treat only the one-dimensional case). ^ is a real rotation angle, and the eigenvalues 
should not depend on it [23 HZI; which is a useful fact for testing convergence. To 
evaluate the nonlinear term IV'P = away from the real coordinate (or x) axis, we 
need to define a generalised complex conjugate which reduces to "^{x) = ■ip{x)* for 
X G M. Applying the complex scaling transformation to 

^ ^ = R{e)i,{x) = e'^l^i){xe'^) , (8) 

we see that ip can be obtained from via the relation: 

if{x) = R{e) (i?(-^)/) * (x) . (9) 
The analytic continuation of Eq. ^ to the complex domain can now be stated as 

H'M< = ' (10) 

with 

H'ii^l] = + ^(^^'') + 9eif,{x)i;l{x) . (11) 

The nonlinear interaction strength is defined here as ge = ge^^^ to compensate for the 
two identical phase factors e'^^^ of if)^ and ■?/' . 

2.2. Numerical solution and propagation algorithm 

In the linear case with (7 = 0, the complex eigenvalue problem of the form Eq. pO|) is 
usually solved by representing the complex-scaled Hamiltonian in a suitable basis and 
final matrix diagonalisation [SH]- For 7^ 0, the corresponding problem to find the 
eigenvalues can be solved only by implicit methods, since H^[^^] explicitly depends on 
the wave function. 

We solved Eq. (fTUj) by searching for the ground-state solution in a self-consistent 
manner. Starting with an initial guess for the wave function ip^^x, t = 0), we evolved in 
real-time the grid representation of ip^{x,t), i.e., 

n 

^%x,t)= 5^ c,(t)x,(x) , (12) 

j=-n 

with the box functions 

, otherwise , 
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and a suitable grid spacing A^.. 

The time-propagation was performed by a sequential application of two different 
integration methods. First, we used a sequence of Crank-Nicholson steps [SO], i-e., 

(1 + iH^At/2) ^/j%x, t + At/2) = (1 - iH^At/2) ^p^x, t - At/2) . (14) 

The Crank-Nicholson method has the advantage of preserving the norm of the wave 
function, but the disadvantage that it treats all modes equally. Since we are interested 
in the ground state, we iterated in a second stage the explicit relation 

/(x, t + At) = (1 - iH%tP^]At) ^%x, t) . (15) 

The latter method, which still corresponds to a real-time integration of the complex 
scaled Gross-Pitaevskii equation, tends to suppress the higher modes [30; and leads 
to a faster stabilisation of the numerical solution of Eq. in comparison with the 
Cranck-Nicholson method ()14|) . For each time step t \—>- t + At, we self-consistently 
solved Eq. (|T5|) by using the left side of Eq. ()15|1 to approximate the nonlinear term 
gei' ip^- Three to five such self-consistent iterations proved sufficient for a stable and 
reliable time propagation. The second derivative appearing in was approximated by 
a finite difference representation (in other words we applied the "forward time centred 
space" representation [SOj to solve the GPE). This leads to a tridiagonal Hamiltonian 
matrix, which significantly simplifies the implementation of both propagators (|T^ and 

For evaluating ip {x,t), we used the method described in detail in which 
produced reliable numerical results also for our Wannier-Stark problem. Briefly 
speaking, we represent ip^{x,t) in a basis set of Gaussians with increasing variance 
for increasing The Gaussian basis is thus well-behaved at the boundaries of our 
grid, which allows us a numerically stable back-rotation to the real domain in x. At 
the end, ip {x, t) is re-expressed again in the grid basis. The necessary matrix- vector 
multiplications are fast since the number of vectors in the Gaussian set can typically be 
chosen much smaller than the number of grid points in the spatial domain. Furthermore, 
the transformation matrices are effectively banded, which reduces the numerical effort 
(note that we computed {x,t) from ip^{x,t) for each time step t t + At to ensure 
stable convergence). 

3. Results and discussion 

In the following, we present our results on the tunnelling rates of resonance states [c.f. 
Eq. (fTUj)] in the Wannier-Stark problem as sketched in Fig. Ufa). Without loss of 
generality we kept fixed the potential depth Vq = 2 in Eq. Q, which corresponds to an 
optical lattice with a maximal amplitude of 16 photon recoil energies [HI EI- We were 
particularly interested in studying the impact of the nonlinear term in Eq. Q on the 
resonant tunnelling peaks of Fig. ^b). Using the method described in the previous 
section we chose 6 = 0.01 . . . 0.02 (where we found stable eigenvalues which are not 
dependent on 6 in this range), and a grid spacing A^ = 0.02 . . . 0.05 for —100 < x < 100. 
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Figure 2. Comparison between the tunnelling rates around the first resonant 
tunnelling peak from Fig. db), obtained for g = by direct diagonalisation of the 
problem (solid line) and by our complex scaling algorithm (circles). 



The integration time step was At = 2.5 X 10"^ for \g\ < 0.2 and reduced to At = 2 x 10~^ 
for larger \g\ > 0.2 and the region of small F < 0.2, while the maximal integration time 
for finding one eigenvalue was tmax = 300. 

As expected, it is very difficult to find the correct eigenvalue close to a resonant 
tunnelling peak because of two reasons: (i) the rates F vary dramatically around the 
peak due to the close degeneracy of the Wannier-Stark levels in adjacent wells, and (ii) 
the rates are rather small 10"^^ < T < 2 x 10~^ at small fields 0.05 < F < 0.2. 

Therefore, we focused on the first (i.e., with m = 1) peak in Fig. ^b), where 
the rates remain F > 10~^. Moreover, we improved the stability of the integration by 
starting with parameters in a stable regime (where we easily found stable, fast converging 
eigenvalues) and adiabatically changing the two parameters F and g into less stable 
parameter regimes. In our case, for Vq = 2, the stable regime is above F > 0.25 for not 
too large nonlinearities \g\ < 0.5 (with optimal stability properties for g = 0). 

We tested our results in three different and independent ways. First, we compared 
them ioT g = with the spectra of the linear Wannier-Stark problem, which can be 
computed by a standard diagonalisation of Hg=o I^H]- Figure El shows the good 
agreement between the data set obtained by our integration of the complex scaled, 
linear version of Eq. ^ and the data presented already in Fig. ^b). 

Secondly, for moderate nonlinearities \g\ < 0.5, we computed for the un-scaled 
problem ^ the survival probability 



with Pc ^ 3 photon recoils (such as to cover the support in momentum space of the initial 




(16) 
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Q 


F 


Ffit 







0.5 


1.94 ± 0.01 X 10"^ 


1.941 X 10"^ 


0.1 


0.5 


2.18 ± 0.01 X 10"^ 


2.180 X 10-^ 





0.25 


7.2 ± 0.1 X lO"'' 


7.2 X 10-^ 


0.1 


0.25 


8.4 ±0.1 X 10"^ 


8.4 X 10-"^ 


0.2 


0.25 


9.7 ±0.1 X 10"^ 


9.7 X 10-^ 


0.25 


0.25 


1.04 ±0.02 X 10-3 


1.04 X 10-3 


0.5 


0.25 


1.45 ±0.03 X 10-3 


1.48 X 10-3 


0.2 


0.15 


3.0 ±0.2 X 10-5 


2.9 X 10-5 


0.2 


0.13125 


5.7 ±0.3 X 10-5 


5.7 X 10-5 



Table 1. Comparison between the tunnelling rates for Vq = 2 obtained by the complex 
scaling method (Fes) and by the integration of the GPE {Tat; integration time up to 
100 Bloch periods; the integration was performed on a large grid that covered the 
full extension of the tunnelled and subsequently accelerated part of the wave function, 
without the use of any cutoff or absorbing boundary conditions). Because of the 
restriction in the integration time, Tfn carries the shown error, whilst the complex 
scaling method allows us to compute the rates Fes with an absolute accuracy of at 
least 10"^ for F > 0.15, and 10"^ for F down to > 0.12. 



state prepared in the spatially periodic lattice potential). Psur(t) was introduced for the 
Wannier-Stark problem in ^H] and characterises the out-coupled loss in momentum 
space, which corresponds to the part of the condensate which has tunnelled through the 
potential. We integrated Eq. Q constantly applying the renormalisation condition 
and computed the survival probability (|T6|l with the wave function 

for the discretised times jt/K {j = 1,2, . . . , K) and large K ^ N. Here %l){x,jt/K) 
denotes the propagated wavefunction immediately before applying the renormalisation 
{ip{x, jt/ K) is renormalised afterwards and propagated up to time (j ± l)t/K). 
■^renorm (p, ^) represents the Fourier transform of the renormalised wavefunction at the 
end of the complete propagation. The decay rates Ffit were obtained by a direct mono- 
exponential fit to the temporal decay of Psur(^)- Table ^ highlights the good agreement 
with the rates computed by the complex-scaling method. 

As a final test of our results, we constantly monitored the quality of the computed 
eigenvalues by evaluating — E)il)^\, which in all cases had to be < 10"*^ for not 

rejecting the eigenvalue. This boundary was chosen such as to be more than two orders 
of magnitude smaller than the smallest tunnelling rates which we computed. 

Our central results are reported now in Fig. IHl There we observe two effects which 
are induced by the presence of the nonlinear interaction term in Eq. (I) the resonant 
tunnelling peak shifts systematically with increasing as a function of the Stark field 
amplitude F. (II) the peak width slightly increases as l^^l increases away from zero. 
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Figure 3. Decay rates F obtained for Va ~ 2 and the range of nonlinearities 
—0.25 < g < 0.25 as a function of the Stark field amphtude F [(a) logarithmic and (b) 
linear scale on y-axis]: g — —0.25 (crosses), g = —0.2 (squares), —0.1 (circles), (solid 
line), 0.1 (pyramids), 0.2 (diamonds; only in (b)), and 0.25 (inverse pyramids; only 
in (b)). The nonlinearity systematically shifts the peak centres and slightly broadens 
the peaks. Their width we defined as the full peak width at half maximum, which is 
marked by the arrows in (b) for the 5 = peak. The dot-dashed line in (a) shows the 
lower bound for converged F > 10~^ on the left side of the resonant tunnelling peaks, 

1. e., for very small F < 0.12, where convergence is very hard to achieve even at small 
5 ~ 0.1 (in this parameter range, i.e., at the left side of the peaks where F changes 
abruptly by about two orders of magnitude, the propagation according to Eqs. H14|l 
and ifT^ becomes unstable). 

The slight broadening goes along with a small increase in the height of the peaks 
with increasing nonlinearity g. Such a destabilisation of the condensate for g > 0, more 
precisely of the decay in the survival probability Psur(^), has already been observed in 
|16j . The ratio of the difference in the height and the difference in the peak width 
(measured at the half of the peak height, see Fig. ^ is roughly constant as a function 
of the nonlinearity < g ^ 0.25. The broadening and the change in height of the 
peak are caused by two different, but simultaneously acting mechanisms. The nonlinear 
mean-field term in Eq. (j2I) partially lifts the degeneracy of the Wannier-Stark levels 
[as sketched in Fig. ^a)] by smearing them out. This qualitatively explains the slight 
broadening of the peak with increasing \g\. In addition, the peak maximum becomes 
systematically larger as g increases in Fig. El because the condensate is destabilised 
(stabilised) by an increasingly repulsive (attractive) nonlinearity. 
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The shift of the peak maximum can be estimated by first-order perturbation theory, 
which predicts the following shift in energy of the levels with respect to the linear case 
with (7 = 0: 



AE^g dx\^g\'\ijg=o\' . (17) 

J -Xc 

For the moderate nonlinearities realised in experiments [HEIIS]! i-e., I*?] ^ 0.5, the 
overlap integral is nearly independent of g, and the major contribution comes from the 
central potential well where the condensate is localised initially. Hence we approximate 
further 

AE^gf dx\ijg=o\\ (18) 

J —TT 

which in turn leads to a shift in the position on the E axis corresponding to AE ^ 271 AE. 
Taking into consideration that the probability in the central well remains normalised 
[condition from Eq. (jS))], Eq. (fTHj) corresponds to the energy shift induced by the 
nonlinearity which was baptised "frequency pulling" in because it leads to a phase 
dispersion of the Bloch oscillations in the wells. With ()18p we arrive at the general 
result that the following ratio is approximately constant: 

—TT^ dx\ipg=o\^0.37 . (19) 

This estimate proved to be valid with a maximal relative deviation of less than 18% 
with respect to the shifts observed in Fig. El Since the integral J^^ is constant 

up to the third digit for all \g\ < 0.25, we conclude that the first-order perturbation 
correction is not enough to describe the shifts more quantitatively. 

While a repulsive mean field {g > 0) enhances the tunnelling rate far away from 
the g = peak, an attractive interaction {g < 0) stabilises the decay sufficiently far 
away from all the peaks. This is the case, e.g., for E > 0.16 in Fig. Efa), where F 
systematically decreases with decreasing g. The symmetric displacement of the peak 
with respect to the sign of g refiects the symmetry of the Bloch band model, in which 
tunnelling from the first excited band back to the ground band can be interpreted as 
the converse process but with a sign change in g JHl- The same qualitative behaviour of 
enhancement {g > 0) and stabilisation {g < 0) was observed in the short-time evolution 
of the three dimensional Wannier-Stark problem JH]- Apart from the conceptional 
difficulty of decomposing a solution of the nonlinear GPE Q into contributions from 
many adjacent wells (see discussion in section E)), the observed washing out of the 
peak structure in Ref. P^Bj is a direct consequence of an effectively moving peak as \g\ 
diminishes monotonously (c.f. Fig. ^ with decreasing density in the wells. 

To conclude this section, we briefly compare our results to other recent works 
which investigated the impact of a mean-fleld interaction of the GPE type on quantum 
mechanical decay processes. We emphasise, however, that such a comparison can only 
be qualitative, since such works [23 I2HI (HHl IM] typically treat much simpler model 
potentials than our Wannier-Stark problem with close level degeneracies. The common 
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feature of our work and the results presented in j2Sl I2H1 ISSl IM| is that an increasing 
nonhnearity typically enhances the decay in the one-dimensional problem. Systematical 
shifts in the chemical potential of resonance states (induced by the interaction term) 
were analytically studied in for a delta-shell potential. Such shifts correspond to 
our perturbative estimate in Eqs. (fTH|) and (fTI?|) . 

Particularly, we compared the results obtained from our method (see section |21) with 
the description of |2H1 where a nonlinear equation was introduced which differs from ours 
[see Eq. flllj) ] in the treatment of the nonlinear term. In the nonhnearity is of the 
form gslip^Y, and we found that such a nonlinear term leads to different decay rates 
than the ones we computed from either our complex scaling method or from the real- 
time integration of the un-scaled GPE ^ and subsequent fits to Psur(^)- We conclude 
that a treatment based on complex scaling - for a condensate within the standard GPE 
description and generally complex-valued wave functions - makes it necessary to use 
the explicit form of as presented in section |21 

There is a growing literature of works on Landau- Zener tunnelling in the presence 
of a mean-field nonhnearity and its impact on the Bloch oscillation problem, see, e.g., 
fSllSSl- In such works a similar systematical stabilisation (for g < 0) ot destabilisation 
(for g > 0) was predicted (and also measured, see [13]) for a single Landau-Zener 
tunnelling event with various approximative models. This corresponds to our results on 
the decay rates which describe directly the initial decay of the condensate via tunnelling, 
i.e., the behaviour of Psur(^) at short times that are not much larger than one Bloch 
period. At and close to the resonant tunnelling peaks the problem is, however, more 
subtle because of the strong interaction of Wannier-Stark levels [see Fig. ^a)], and such 
a case was not treated in [l31 EHj ■ 

4. Conclusion 

To summarise, we presented a method to numerically compute precise decay rates for 
tunnelling problems within the framework of the Gross-Pitaevskii equation. We adapted 
and improved the technique developed by one us in for the more complicated scenario 
of resonant tunnelling in the Wannier-Stark problem. We showed that the mean- 
field nonhnearity leads to experimentally observable modifications in the tunnelling of 
resonance states from the periodical potential wells, even in a regime where the kinetic 
and the periodic potential terms still dominate the dynamics. The broadening and the 
shift of the resonant tunnelling peaks define clear signatures for nonhnearity induced 
effects. 

Our method can be extended - with further system-specific improvements in 
the propagation algorithm - to treat even more complicated problems appearing in 
experiments with Rose condensates, e.g., the transport of coherent matter within atomic 
wave guides. Finally, we can readily extend the proposed method to three spatial 
dimensions, with the only drawback of much larger numerical effort. 
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